Univariate Arima Models

Predicting COVID-19 Deaths (STAT 390)

Author

Ryan Nguyen

1. Load Libraries, Seed, & Data

Code
# Load necessary libraries
library(forecast)
library(tseries)
library(tidyverse)
library(tidymodels)
library(kableExtra)
library(modeltime)
library(timetk)
library(cowplot)

# setting the seed
set.seed(123)

# Loading Data  ----------------------------------------------------------------
daily_deaths <- read_csv("data/new_daily_deaths.csv") %>% 
  select(daily_deaths)

east <- read_csv("data/east_daily.csv")

midwest <- read_csv("data/midwest_daily.csv")

south <- read_csv("data/south_daily.csv")

west <- read_csv("data/west_daily.csv")

2. Data Splitting

Code
# Data Splitting ---------------------------------------------------------------
east_splits <- time_series_split(east, initial = "1022 days", assess = "114 days")
east_train <- training(east_splits)
east_test <- testing(east_splits)

midwest_splits <- time_series_split(midwest, initial = "1022 days", assess = "114 days")
midwest_train <- training(midwest_splits)
midwest_test <- testing(midwest_splits)

south_splits <- time_series_split(south, initial = "1022 days", assess = "114 days")
south_train <- training(south_splits)
south_test <- testing(south_splits)

west_splits <- time_series_split(west, initial = "1022 days", assess = "114 days")
west_train <- training(west_splits)
west_test <- testing(west_splits)

## Convert Into Time Series ----------------------------------------------------
east_train_ts <- ts(east_train %>%
                      select(daily_deaths))
east_test_ts <- ts(east_test %>%
                     select(daily_deaths))
east_ts <- ts(east %>% 
                select(daily_deaths))

midwest_train_ts <- ts(midwest_train %>%
                      select(daily_deaths))
midwest_test_ts <- ts(midwest_test %>%
                     select(daily_deaths))
midwest_ts <- ts(midwest %>% 
                select(daily_deaths))

south_train_ts <- ts(south_train %>%
                    select(daily_deaths))
south_test_ts <- ts(south_test %>%
                   select(daily_deaths))
south_ts <- ts(south %>%
                 select(daily_deaths))

west_train_ts <- ts(west_train %>%
                   select(daily_deaths))
west_test_ts <- ts(west_test %>%
                  select(daily_deaths))
west_ts <- ts(west %>%
                select(daily_deaths))

## Plotting Training and Testing -----------------------------------------------
par(mfrow=c(2,2))
ggplot() +
  geom_line(data = east_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = east_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = east_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = east_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "East Region Training vs. Testing Data")

Code
ggplot() +
  geom_line(data = midwest_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = midwest_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = midwest_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = midwest_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "Midwest Region Training vs. Testing Data")

Code
ggplot() +
  geom_line(data = south_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = south_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = south_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = south_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "South Region Training vs. Testing Data")

Code
ggplot() +
  geom_line(data = west_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = west_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = west_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = west_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "West Region Training vs. Testing Data")

3. Testing for Stationarity

Code
### Using Augmented Dickey-Fuller Test -----------------------------------------
adf_east <- adf.test(east_ts)
print(adf_east) # p-value is 0.01 < 0.05, implying it is stationary

    Augmented Dickey-Fuller Test

data:  east_ts
Dickey-Fuller = -4.5566, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_midwest <- adf.test(midwest_ts)
print(adf_midwest) # p-value is 0.325 > 0.05, implying it is NOT stationary

    Augmented Dickey-Fuller Test

data:  midwest_ts
Dickey-Fuller = -2.5984, Lag order = 10, p-value = 0.325
alternative hypothesis: stationary
Code
midwest_ts <- diff(midwest_ts)
adf_midwest <- adf.test(midwest_ts)
print(adf_midwest) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff

    Augmented Dickey-Fuller Test

data:  midwest_ts
Dickey-Fuller = -13.43, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_south <- adf.test(south_ts)
print(adf_south) # p-value is 0.3018 > 0.05, implying it is NOT stationary

    Augmented Dickey-Fuller Test

data:  south_ts
Dickey-Fuller = -2.6533, Lag order = 10, p-value = 0.3018
alternative hypothesis: stationary
Code
south_ts <- diff(south_ts)
adf_south <- adf.test(south_ts)
print(adf_south) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff

    Augmented Dickey-Fuller Test

data:  south_ts
Dickey-Fuller = -11.36, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_west <- adf.test(west_ts)
print(adf_west) # p-value is 0.4098 > 0.05, implying it is NOT stationary

    Augmented Dickey-Fuller Test

data:  west_ts
Dickey-Fuller = -2.3981, Lag order = 10, p-value = 0.4098
alternative hypothesis: stationary
Code
west_ts <- diff(west_ts)
adf_west <- adf.test(west_ts)
print(adf_west)  # p-value is 0.01 < 0.05, implying it is stationary with 1 diff

    Augmented Dickey-Fuller Test

data:  west_ts
Dickey-Fuller = -12.247, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
### Making Testing/Training Stationary -----------------------------------------
midwest_train_ts <- diff(midwest_train_ts)
midwest_test_ts <- diff(midwest_test_ts)

midwest_train_ts <- as.data.frame(midwest_train_ts) %>% 
  mutate(daily_deaths = case_when((daily_deaths < 0) ~ 0,
                                  .default = daily_deaths))
midwest_test_ts <- as.data.frame(midwest_test_ts) %>% 
  mutate(daily_deaths = case_when((daily_deaths < 0) ~ 0,
                                  .default = daily_deaths))

midwest_train_ts <- ts(midwest_train_ts)
midwest_test_ts <- ts(midwest_test_ts)

south_train_ts <- diff(south_train_ts)
south_test_ts <- diff(south_test_ts)

south_train_ts <- as.data.frame(south_train_ts) %>% 
  mutate(daily_deaths = case_when((daily_deaths < 0) ~ 0,
                                  .default = daily_deaths))
south_test_ts <- as.data.frame(south_test_ts) %>% 
  mutate(daily_deaths = case_when((daily_deaths < 0) ~ 0,
                                  .default = daily_deaths))

south_train_ts <- ts(south_train_ts)
south_test_ts <- ts(south_test_ts)

west_train_ts <- diff(west_train_ts)
west_test_ts <- diff(west_test_ts)

west_train_ts <- as.data.frame(west_train_ts) %>% 
  mutate(daily_deaths = case_when((daily_deaths < 0) ~ 0,
                                  .default = daily_deaths))
west_test_ts <- as.data.frame(west_test_ts) %>% 
  mutate(daily_deaths = case_when((daily_deaths < 0) ~ 0,
                                  .default = daily_deaths))

west_train_ts <- ts(west_train_ts)
west_test_ts <- ts(west_test_ts)

4. ACF & PACF Plots

Code
## ACF & PACF Plots ------------------------------------------------------------
par(mfrow=c(3, 1))
acf(east_ts, main = "East")
pacf(east_ts, main = "East")
plot(east_ts, main = "Daily Deaths: East")

Code
acf(midwest_ts, main = "Midwest")
pacf(midwest_ts, main = "Midwest")
plot(midwest_ts, main = "Daily Deaths: Midwest")

Code
acf(south_ts, main = "South")
pacf(south_ts, main = "South")
plot(south_ts, main = "Daily Deaths: South")

Code
acf(west_ts, main = "West")
pacf(west_ts, main = "West")
plot(south_ts, main = "Daily Deaths: South")

Code
plot(west_ts, main = "Daily Deaths: West")

5. Time Series Decomposition

Code
# total daily death decomposition
daily_ts <- ts(daily_deaths, frequency = 365)

decomposed_ts <- decompose(daily_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(daily_ts, main = "Total Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_ts$trend, main = "Trend Component", ylab = "Value", col = "red")

Code
plot(decomposed_ts$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_ts$random, main = "Residual Component", ylab = "Value", col = "purple")

Code
# East Decomposition -----------------------------------------------------------
east_daily <- east %>% 
  select(daily_deaths)

east_ts <- ts(east_daily, frequency = 365)

decomposed_east <- decompose(east_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(east_ts, main = "East Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_east$trend, main = "Trend Component", ylab = "Value", col = "red")

Code
plot(decomposed_east$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_east$random, main = "Residual Component", ylab = "Value", col = "purple")

Code
## Midwest Decomposition ----------------------------------------------------------
midwest_daily <- midwest %>% 
  select(daily_deaths)

midwest_ts <- ts(midwest_daily, frequency = 365)

decomposed_midwest <- decompose(midwest_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(midwest_ts, main = "Midwest Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_midwest$trend, main = "Trend Component", ylab = "Value", col = "red")

Code
plot(decomposed_midwest$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_midwest$random, main = "Residual Component", ylab = "Value", col = "purple")

Code
## South Decomposition ----------------------------------------------------------
south_daily <- south %>% 
  select(daily_deaths)

south_ts <- ts(south_daily, frequency = 365)

decomposed_south <- decompose(south_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(south_ts, main = "South Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_south$trend, main = "Trend Component", ylab = "Value", col = "red")

Code
plot(decomposed_south$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_south$random, main = "Residual Component", ylab = "Value", col = "purple")

Code
## West Decomposition ----------------------------------------------------------
west_daily <- west %>% 
  select(daily_deaths)

west_ts <- ts(west_daily, frequency = 365)

decomposed_west <- decompose(west_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(west_ts, main = "West Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_west$trend, main = "Trend Component", ylab = "Value", col = "red")

Code
plot(decomposed_west$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_west$random, main = "Residual Component", ylab = "Value", col = "purple")

6. ARIMA Model

6.1 Initial Model

Code
# Building ARIMA Model ---------------------------------------------------------
east_arima <- arima(east_train_ts, order=c(1,0,1))
midwest_arima <- arima(midwest_train_ts, order=c(1,1,1))
south_arima <- arima(south_train_ts, order=c(1,1,1))
west_arima <- arima(west_train_ts, order=c(1,1,1))

## Summary & Forecast ---------------------------------------------------------
summary(east_arima)

Call:
arima(x = east_train_ts, order = c(1, 0, 1))

Coefficients:
         ar1      ma1  intercept
      0.9176  -0.0388   197.3738
s.e.  0.0147   0.0493    39.7998

sigma^2 estimated as 12153:  log likelihood = -6257.17,  aic = 12522.33

Training set error measures:
                    ME     RMSE      MAE  MPE MAPE      MASE        ACF1
Training set 0.1151898 110.2407 64.66739 -Inf  Inf 0.9918341 0.008515868
Code
summary(midwest_arima)

Call:
arima(x = midwest_train_ts, order = c(1, 1, 1))

Coefficients:
          ar1      ma1
      -0.2237  -0.9343
s.e.   0.0313   0.0090

sigma^2 estimated as 16650:  log likelihood = -6405.86,  aic = 12817.73

Training set error measures:
                    ME     RMSE     MAE MPE MAPE      MASE        ACF1
Training set 0.7797006 128.9735 78.3255 NaN  Inf 0.7410171 -0.06964661
Code
summary(south_arima)

Call:
arima(x = south_train_ts, order = c(1, 1, 1))

Coefficients:
          ar1      ma1
      -0.1913  -0.9453
s.e.   0.0317   0.0095

sigma^2 estimated as 44608:  log likelihood = -6908.51,  aic = 13823.02

Training set error measures:
                   ME     RMSE     MAE MPE MAPE      MASE        ACF1
Training set 1.277928 211.1023 129.556 NaN  Inf 0.7498262 -0.01884038
Code
summary(west_arima)

Call:
arima(x = west_train_ts, order = c(1, 1, 1))

Coefficients:
          ar1      ma1
      -0.0842  -0.9603
s.e.   0.0319   0.0072

sigma^2 estimated as 10324:  log likelihood = -6162.22,  aic = 12330.44

Training set error measures:
                    ME    RMSE      MAE  MPE MAPE      MASE         ACF1
Training set 0.6807731 101.558 62.19197 -Inf  Inf 0.7886397 -0.009430079
Code
east_values <- forecast(east_arima, h = length(east_test_ts))
midwest_values <- forecast(midwest_arima, h = length(midwest_test_ts))
south_values <- forecast(south_arima, h = length(south_test_ts))
west_values <- forecast(west_arima, h = length(west_test_ts))

6.2 Accuracy

Code
## Evaluating Accuracy ---------------------------------------------------------
east_initial_arima_acc <- forecast::accuracy(east_values)

east_initial_arima_acc %>% 
  kbl(caption = "<center>Initial East Accuracy<center>") %>% 
  kable_styling()
Initial East Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1151898 110.2407 64.66739 -Inf Inf 0.9918341 0.0085159
Code
### East Region
east_test_plot <- as.data.frame(east_test_ts) %>% 
  mutate(num = seq(1023, (1023+113)))

east_initial_arima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Initial ARIMA",
       y = "Daily Deaths")

east_initial_arima_plot

Code
### Midwest Region

midwest_initial_arima_acc <- forecast::accuracy(midwest_values)

midwest_initial_arima_acc %>% 
  kbl(caption = "<center>Initial Midwest Accuracy<center>") %>% 
  kable_styling()
Initial Midwest Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.7797006 128.9735 78.3255 NaN Inf 0.7410171 -0.0696466
Code
midwest_test_plot <- as.data.frame(midwest_test_ts) %>% 
  mutate(num = seq(1023, (1023+112)))

midwest_initial_arima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Initial ARIMA",
       y = "Daily Deaths")

midwest_initial_arima_plot

Code
### South Region

south_initial_arima_acc <- forecast::accuracy(south_values)

south_initial_arima_acc %>% 
  kbl(caption = "<center>Initial South Accuracy<center>") %>% 
  kable_styling()
Initial South Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 1.277928 211.1023 129.556 NaN Inf 0.7498262 -0.0188404
Code
south_test_plot <- as.data.frame(south_test_ts) %>% 
  mutate(num = seq(1023, (1023+112)))

south_initial_arima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Initial ARIMA",
       y = "Daily Deaths")

south_initial_arima_plot 

Code
### West Region

west_initial_arima_acc <- forecast::accuracy(west_values)

west_initial_arima_acc %>% 
  kbl(caption = "<center>Initial West Accuracy<center>") %>% 
  kable_styling()
Initial West Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.6807731 101.558 62.19197 -Inf Inf 0.7886397 -0.0094301
Code
west_test_plot <- as.data.frame(west_test_ts) %>% 
  mutate(num = seq(1023, (1023+112)))

west_initial_arima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Initial ARIMA",
       y = "Daily Deaths") 

west_initial_arima_plot

6.4 Rebuilding ARIMA Model

Code
# Rebuilding ARIMA Models ------------------------------------------------------
east_arima <- arima(east_train_ts, order=c(5,0,5))
midwest_arima <- arima(midwest_train_ts, order=c(5,1,5))
south_arima <- arima(south_train_ts, order=c(5,1,5))
west_arima <- arima(west_train_ts, order=c(5,1,5), method = "ML")


east_values <- forecast(east_arima, h = length(east_test_ts))
midwest_values <- forecast(midwest_arima, h = length(midwest_test_ts))
south_values <- forecast(south_arima, h = length(south_test_ts))
west_values <- forecast(west_arima, h = length(west_test_ts))

## Evaluating Accuracy & Plots -------------------------------------------------

east_tuned_arima_acc <- forecast::accuracy(east_values)

east_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned East Accuracy<center>") %>% 
  kable_styling()
Tuned East Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1095068 80.35755 47.18823 NaN Inf 0.723748 -0.0254917
Code
### East Region
east_tuned_arima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Tuned ARIMA",
       y = "Daily Deaths")

east_tuned_arima_plot

Code
### Midwest Region

midwest_tuned_arima_acc <- forecast::accuracy(midwest_values)

midwest_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned Midwest Accuracy<center>") %>% 
  kable_styling()
Tuned Midwest Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.4882952 96.70073 54.43655 NaN Inf 0.5150099 -0.0118025
Code
midwest_tuned_arima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Tuned ARIMA",
       y = "Daily Deaths")

midwest_tuned_arima_plot

Code
### South Region

south_tuned_arima_acc <- forecast::accuracy(south_values)

south_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned South Accuracy<center>") %>% 
  kable_styling()
Tuned South Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.7086084 172.0968 94.70831 NaN Inf 0.5481396 -0.0036251
Code
south_tuned_arima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Tuned ARIMA",
       y = "Daily Deaths")

south_tuned_arima_plot 

Code
### West Region
west_tuned_arima_acc <- forecast::accuracy(west_values)

west_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned West Accuracy<center>") %>% 
  kable_styling()
Tuned West Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.4533154 79.29443 46.88033 NaN Inf 0.5944769 -0.0091414
Code
west_tuned_arima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Tuned ARIMA",
       y = "Daily Deaths") 

west_tuned_arima_plot

Code
# saving results
save(east_tuned_arima_acc, east_initial_arima_acc, east_initial_arima_plot, east_tuned_arima_plot,
     midwest_tuned_arima_acc, midwest_initial_arima_acc, midwest_initial_arima_plot, midwest_tuned_arima_plot,
     south_tuned_arima_acc, south_initial_arima_acc, south_initial_arima_plot, south_tuned_arima_plot,
     west_tuned_arima_acc, west_initial_arima_acc, west_initial_arima_plot, west_tuned_arima_plot,
     file = "results/uni_arima.rda"
     )

7. SARIMA Model

7.1 Initial SARIMA Model

Code
# Building SARIMA model  ------------------------------------------------------
east_sarima <- arima(east_train_ts, seasonal = list(order = c(1,1,1)))
midwest_sarima <- arima(midwest_train_ts, seasonal = list(order = c(1,1,1)))
south_sarima <- arima(south_train_ts, seasonal = list(order = c(1,1,1)))
west_sarima <- arima(west_train_ts, seasonal = list(order = c(1,1,1)))

# Summary of the SARIMA model
summary(east_sarima)

Call:
arima(x = east_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
        sar1     sma1
      0.4706  -0.8016
s.e.  0.0382   0.0211

sigma^2 estimated as 11156:  log likelihood = -6206.64,  aic = 12419.27

Training set error measures:
                    ME     RMSE     MAE  MPE MAPE      MASE      ACF1
Training set 0.1606608 105.5715 61.2767 -Inf  Inf 0.9398295 0.1140208
Code
summary(midwest_sarima)

Call:
arima(x = midwest_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
         sar1     sma1
      -0.2237  -0.9343
s.e.   0.0313   0.0090

sigma^2 estimated as 16650:  log likelihood = -6405.86,  aic = 12817.73

Training set error measures:
                    ME     RMSE     MAE MPE MAPE      MASE        ACF1
Training set 0.7797006 128.9735 78.3255 NaN  Inf 0.7410171 -0.06964661
Code
summary(south_sarima)

Call:
arima(x = south_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
         sar1     sma1
      -0.1913  -0.9453
s.e.   0.0317   0.0095

sigma^2 estimated as 44608:  log likelihood = -6908.51,  aic = 13823.02

Training set error measures:
                   ME     RMSE     MAE MPE MAPE      MASE        ACF1
Training set 1.277928 211.1023 129.556 NaN  Inf 0.7498262 -0.01884038
Code
summary(west_sarima)

Call:
arima(x = west_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
         sar1     sma1
      -0.0842  -0.9603
s.e.   0.0319   0.0072

sigma^2 estimated as 10324:  log likelihood = -6162.22,  aic = 12330.44

Training set error measures:
                    ME    RMSE      MAE  MPE MAPE      MASE         ACF1
Training set 0.6807731 101.558 62.19197 -Inf  Inf 0.7886397 -0.009430079
Code
## Forecasting with the SARIMA model -------------------------------------------
east_values <- forecast(east_sarima, h = length(east_test_ts))
midwest_values <- forecast(midwest_sarima, h = length(midwest_test_ts))
south_values <- forecast(south_sarima, h = length(south_test_ts))
west_values <- forecast(west_sarima, h = length(west_test_ts))

7.2 Accuracy

Code
## Plotting & Accuracy

## East Region
east_initial_sarima_acc <- forecast::accuracy(east_values)

east_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA East Accuracy<center>") %>% 
  kable_styling()
Initial SARIMA East Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1606608 105.5715 61.2767 -Inf Inf 0.9398295 0.1140208
Code
east_initial_sarima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Initial SARIMA",
       y = "Daily Deaths")

east_initial_sarima_plot

Code
## Midwest Region

midwest_initial_sarima_acc <- forecast::accuracy(midwest_values)

midwest_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA Midwest Accuracy<center>") %>% 
  kable_styling()
Initial SARIMA Midwest Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.7797006 128.9735 78.3255 NaN Inf 0.7410171 -0.0696466
Code
midwest_initial_sarima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Initial SARIMA",
       y = "Daily Deaths")

midwest_initial_sarima_plot

Code
## South Region
south_initial_sarima_acc <- forecast::accuracy(south_values)

south_initial_sarima_acc %>%
  kbl(caption = "<center>Initial SARIMA South Accuracy<center>") %>%
  kable_styling()
Initial SARIMA South Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 1.277928 211.1023 129.556 NaN Inf 0.7498262 -0.0188404
Code
south_initial_sarima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Initial SARIMA",
       y = "Daily Deaths")

south_initial_sarima_plot

Code
## West Region

west_initial_sarima_acc <- forecast::accuracy(west_values)

west_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA West Accuracy<center>") %>% 
  kable_styling()
Initial SARIMA West Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.6807731 101.558 62.19197 -Inf Inf 0.7886397 -0.0094301
Code
west_initial_sarima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Initial SARIMA",
       y = "Daily Deaths")

west_initial_sarima_plot

7.3 Grid Search

Code
## Define the parameter grids --------------------------------------------------
p_grid <- 0:5  # AR parameter
d_grid <- 0:0  # Differencing parameter
q_grid <- 0:5  # MA parameter
P_grid <- 0:1  # Seasonal AR parameter
D_grid <- 0:1  # Seasonal differencing parameter
Q_grid <- 0:1  # Seasonal MA parameter
s <- 12        # Seasonal period



## East Grid Search ------------------------------------------------------------

# Create East Region
east_train <- east_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
east_results <- data.frame(order = character(),
                      seasonal = character(),
                      AIC = numeric(),
                      stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(east_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                east_results <- rbind(east_results, data.frame(order = paste(order, collapse = ", "),
                                                     seasonal = paste(seasonal, collapse = ", "),
                                                     AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}

d_grid <- 0:1  # Differencing parameter

## Midwest Grid Search ------------------------------------------------------------

# Create Midwest Region
midwest_train <- midwest_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
midwest_results <- data.frame(order = character(),
                           seasonal = character(),
                           AIC = numeric(),
                           stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(midwest_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                midwest_results <- rbind(midwest_results, data.frame(order = paste(order, collapse = ", "),
                                                               seasonal = paste(seasonal, collapse = ", "),
                                                               AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}

## South Grid Search ------------------------------------------------------------

# Create South Region
south_train <- south_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
south_results <- data.frame(order = character(),
                           seasonal = character(),
                           AIC = numeric(),
                           stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(south_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                south_results <- rbind(south_results, data.frame(order = paste(order, collapse = ", "),
                                                               seasonal = paste(seasonal, collapse = ", "),
                                                               AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}


## West Grid Search ------------------------------------------------------------

# Create West Region
west_train <- west_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
west_results <- data.frame(order = character(),
                           seasonal = character(),
                           AIC = numeric(),
                           stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(west_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                west_results <- rbind(west_results, data.frame(order = paste(order, collapse = ", "),
                                                               seasonal = paste(seasonal, collapse = ", "),
                                                               AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}
Code
save(east_results, midwest_results, south_results, west_results, file = "results/sarima_grid.rda")
Code
load("results/sarima_grid.rda")

Best East SARIMA

Code
# Find the best model (minimum AIC)
best_east_model <- east_results[which.min(east_results$AIC), ]
print(best_east_model)
      order    seasonal      AIC
213 5, 0, 4 0, 1, 1, 12 11865.53

Best Midwest SARIMA

Code
# Midwest Results
# Find the best model (minimum AIC)
best_midwest_model <- midwest_results[which.min(midwest_results$AIC), ]
print(best_midwest_model)
      order    seasonal      AIC
385 5, 1, 5 0, 0, 1, 12 12819.82

Best South SARIMA

Code
# South Results
# Find the best model (minimum AIC)
best_south_model <- south_results[which.min(south_results$AIC), ]
print(best_south_model)
      order    seasonal      AIC
402 5, 1, 5 0, 0, 0, 12 13698.93

Best West SARIMA

Code
# Find the best model (minimum AIC)
best_west_model <- west_results[which.min(west_results$AIC), ]
print(best_west_model)
      order    seasonal      AIC
385 5, 0, 5 1, 1, 0, 12 12141.63

7.4 Rebuilding SARIMA

Code
# Rebuilding SARIMA Models -----------------------------------------------------
east_sarima <- arima(east_train_ts, order = c(5, 0, 4), 
                     seasonal = list(order = c(0, 1, 1), period = 12))
midwest_sarima <- arima(midwest_train_ts, order = c(5, 1, 5),
                        seasonal = list(order = c(0, 0, 1), period = 12))
south_sarima <- arima(south_train_ts, order = c(5, 1, 5),
                           seasonal = list(order = c(0, 0, 0), period = 12))
west_sarima <- arima(west_train_ts, order = c(5, 0, 5),
                     seasonal = list(order = c(1, 1, 0), period = 12))


## Forecasting with the NEW SARIMA model ---------------------------------------
east_values <- forecast(east_sarima, h = length(east_test_ts))
midwest_values <- forecast(midwest_sarima, h = length(midwest_test_ts))
south_values <- forecast(south_sarima, h = length(south_test_ts))
west_values <- forecast(west_sarima, h = length(west_test_ts))
Code
## Plotting & Accuracy

## East Region
east_tuned_sarima_acc <- forecast::accuracy(east_values)

east_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA East Accuracy<center>") %>% 
  kable_styling()
Tuned SARIMA East Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1031443 87.81972 53.94287 NaN Inf 0.8273471 -0.0998164
Code
east_tuned_sarima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Tuned SARIMA",
       y = "Daily Deaths")

east_tuned_sarima_plot

Code
## Midwest Region
midwest_tuned_sarima_acc <- forecast::accuracy(midwest_values)

midwest_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA Midwest Accuracy<center>") %>% 
  kable_styling()
Tuned SARIMA Midwest Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.4955934 106.0021 63.95655 NaN Inf 0.6050762 0.0305554
Code
midwest_tuned_sarima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Tuned SARIMA",
       y = "Daily Deaths")

midwest_tuned_sarima_plot

Code
## South Region
south_tuned_sarima_acc <- forecast::accuracy(south_values)

south_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA South Accuracy<center>") %>% 
  kable_styling()
Tuned SARIMA South Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.7086084 172.0968 94.70831 NaN Inf 0.5481396 -0.0036251
Code
south_tuned_sarima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Tuned SARIMA",
       y = "Daily Deaths")

south_tuned_sarima_plot

Code
## West Region
west_tuned_sarima_acc <- forecast::accuracy(west_values)

west_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA West Accuracy<center>") %>% 
  kable_styling()
Tuned SARIMA West Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.681019 89.05196 55.78604 NaN Inf 0.7074078 0.0029928
Code
west_tuned_sarima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Tuned SARIMA",
       y = "Daily Deaths")

west_tuned_sarima_plot

Code
# saving results
save(east_tuned_sarima_acc, east_initial_sarima_acc, east_initial_sarima_plot, east_tuned_sarima_plot,
     midwest_tuned_sarima_acc, midwest_initial_sarima_acc, midwest_initial_sarima_plot, midwest_tuned_sarima_plot,
     south_tuned_sarima_acc, south_initial_sarima_acc, south_initial_sarima_plot, south_tuned_sarima_plot,
     west_tuned_sarima_acc, west_initial_sarima_acc, west_initial_sarima_plot, west_tuned_sarima_plot,
     file = "results/uni_sarima.rda"
     )

8. Auto ARIMA Model

8.1 Building AutoARIMA

Code
# Building AutoARIMA Model -----------------------------------------------------
east_arima <- auto.arima(east_train_ts, seasonal = TRUE)
midwest_arima <- auto.arima(midwest_train_ts, seasonal = TRUE)
south_arima <- auto.arima(south_train_ts, seasonal = TRUE)
west_arima <- auto.arima(west_train_ts, seasonal = TRUE)

# Summary of the ARIMA model
summary(east_arima)
Series: east_train_ts 
ARIMA(5,1,2) 

Coefficients:
         ar1      ar2      ar3      ar4      ar5      ma1     ma2
      0.5142  -0.7304  -0.1394  -0.3479  -0.2207  -1.0295  0.7353
s.e.  0.0480   0.0344   0.0426   0.0329   0.0418   0.0410  0.0265

sigma^2 = 6861:  log likelihood = -5956.39
AIC=11928.79   AICc=11928.93   BIC=11968.22

Training set error measures:
                    ME     RMSE     MAE MPE MAPE      MASE       ACF1
Training set 0.1506172 82.50546 47.1043 NaN  Inf 0.7224607 0.01707355
Code
summary(midwest_arima)
Series: midwest_train_ts 
ARIMA(5,1,4) 

Coefficients:
          ar1      ar2      ar3      ar4      ar5      ma1     ma2      ma3
      -0.0349  -0.7605  -0.0404  -0.1320  -0.4194  -1.1696  0.9929  -0.9722
s.e.   0.0461   0.0370   0.0572   0.0371   0.0377   0.0441  0.0492   0.0459
         ma4
      0.4244
s.e.  0.0381

sigma^2 = 10525:  log likelihood = -6168.66
AIC=12357.32   AICc=12357.54   BIC=12406.6

Training set error measures:
                   ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set 0.385353 102.0896 59.56282 NaN  Inf 0.5635082 -0.05532994
Code
summary(south_arima)
Series: south_train_ts 
ARIMA(1,1,2) 

Coefficients:
         ar1      ma1     ma2
      0.4799  -1.6879  0.7236
s.e.  0.0578   0.0428  0.0405

sigma^2 = 42636:  log likelihood = -6884.08
AIC=13776.16   AICc=13776.2   BIC=13795.87

Training set error measures:
                    ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set 0.7155848 206.0793 126.7407 NaN  Inf 0.7335323 0.008694211
Code
summary(west_arima)
Series: west_train_ts 
ARIMA(5,1,3) 

Coefficients:
         ar1      ar2      ar3      ar4      ar5      ma1     ma2      ma3
      0.4572  -0.6143  -0.2778  -0.2654  -0.2584  -1.8532  1.6316  -0.6651
s.e.  0.0588   0.0355   0.0387   0.0340   0.0435   0.0544  0.0868   0.0458

sigma^2 = 6300:  log likelihood = -5907.95
AIC=11833.91   AICc=11834.09   BIC=11878.26

Training set error measures:
                  ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set 0.44309 79.02328 45.78169 NaN  Inf 0.5805454 0.003263827
Code
## Forecasting Test Data -------------------------------------------------------
east_values <- forecast(east_arima, h = length(east_test_ts))
midwest_values <- forecast(midwest_arima, h = length(midwest_test_ts))
south_values <- forecast(south_arima, h = length(south_test_ts))
west_values <- forecast(west_arima, h = length(west_test_ts))

8.2 Accuracy

Code
## Plotting & Accuracy

## East Region
east_auto_acc <- forecast::accuracy(east_values)

east_auto_acc  %>% 
  kbl(caption = "<center>Auto ARIMA East Accuracy<center>") %>% 
  kable_styling()
Auto ARIMA East Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1506172 82.50546 47.1043 NaN Inf 0.7224607 0.0170736
Code
east_auto_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Auto ARIMA",
       y = "Daily Deaths")

east_auto_plot 

Code
## Midwest Region
midwest_auto_acc <- forecast::accuracy(midwest_values)

midwest_auto_acc %>% 
  kbl(caption = "<center>Auto ARIMA Midwest Accuracy<center>") %>% 
  kable_styling()
Auto ARIMA Midwest Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.385353 102.0896 59.56282 NaN Inf 0.5635082 -0.0553299
Code
midwest_auto_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Auto ARIMA",
       y = "Daily Deaths")

midwest_auto_plot 

Code
## South Region
south_auto_acc <- forecast::accuracy(south_values)

south_auto_acc %>% 
  kbl(caption = "<center>Auto ARIMA South Accuracy<center>") %>% 
  kable_styling()
Auto ARIMA South Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.7155848 206.0793 126.7407 NaN Inf 0.7335323 0.0086942
Code
south_auto_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Auto ARIMA",
       y = "Daily Deaths")

south_auto_plot

Code
## West Region
west_auto_acc <- forecast::accuracy(west_values)

west_auto_acc %>% 
  kbl(caption = "<center>Auto ARIMA West Accuracy<center>") %>% 
  kable_styling()
Auto ARIMA West Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.44309 79.02328 45.78169 NaN Inf 0.5805454 0.0032638
Code
west_auto_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Auto ARIMA",
       y = "Daily Deaths")

west_auto_plot

Code
# saving results
save(east_auto_plot,
     midwest_auto_plot,
     south_auto_plot,
     west_auto_plot,
     east_auto_acc,
     midwest_auto_acc,
     south_auto_acc,
     west_auto_acc,
     file = "results/uni_autoarima.rda")

9. Compare Arima Models

Code
load("results/uni_arima.rda")
load("results/uni_sarima.rda")
load("results/uni_autoarima.rda")

East Results

Code
east1 <- east_initial_arima_acc %>% 
  as.data.frame()

east2 <- east_tuned_arima_acc %>% 
  as.data.frame()

east3 <- east_initial_sarima_acc%>% 
  as.data.frame() 

east4 <- east_tuned_sarima_acc%>% 
  as.data.frame()

east5 <- east_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", east1$MAE, east1$MASE,
  "tuned arima", east2$MAE, east2$MASE,
  "initial sarima", east3$MAE, east3$MASE,
  "tuned sarima", east4$MAE, east4$MASE,
  "auto arima", east5$MAE, east5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()
model mae mase
initial arima 64.66739 0.9918341
tuned arima 47.18823 0.7237480
initial sarima 61.27670 0.9398295
tuned sarima 53.94287 0.8273471
auto arima 47.10430 0.7224607

Midwest Results

Code
midwest1 <- midwest_initial_arima_acc %>% 
  as.data.frame()

midwest2 <- midwest_tuned_arima_acc %>% 
  as.data.frame()

midwest3 <- midwest_initial_sarima_acc%>% 
  as.data.frame() 

midwest4 <- midwest_tuned_sarima_acc%>% 
  as.data.frame()

midwest5 <- midwest_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", midwest1$MAE, midwest1$MASE,
  "tuned arima", midwest2$MAE, midwest2$MASE,
  "initial sarima", midwest3$MAE, midwest3$MASE,
  "tuned sarima", midwest4$MAE, midwest4$MASE,
  "auto arima", midwest5$MAE, midwest5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()
model mae mase
initial arima 78.32550 0.7410171
tuned arima 54.43655 0.5150099
initial sarima 78.32550 0.7410171
tuned sarima 63.95655 0.6050762
auto arima 59.56282 0.5635082

South Results

Code
south1 <- south_initial_arima_acc %>% 
  as.data.frame()

south2 <- south_tuned_arima_acc %>% 
  as.data.frame()

south3 <- south_initial_sarima_acc%>% 
  as.data.frame() 

south4 <- south_tuned_sarima_acc%>% 
  as.data.frame()

south5 <- south_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", south1$MAE, south1$MASE,
  "tuned arima", south2$MAE, south2$MASE,
  "initial sarima", south3$MAE, south3$MASE,
  "tuned sarima", south4$MAE, south4$MASE,
  "auto arima", south5$MAE, south5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()
model mae mase
initial arima 129.55599 0.7498262
tuned arima 94.70831 0.5481396
initial sarima 129.55599 0.7498262
tuned sarima 94.70831 0.5481396
auto arima 126.74072 0.7335323

West Results

Code
west1 <- west_initial_arima_acc %>% 
  as.data.frame()

west2 <- west_tuned_arima_acc %>% 
  as.data.frame()

west3 <- west_initial_sarima_acc%>% 
  as.data.frame() 

west4 <- west_tuned_sarima_acc%>% 
  as.data.frame()

west5 <- west_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", west1$MAE, west1$MASE,
  "tuned arima", west2$MAE, west2$MASE,
  "initial sarima", west3$MAE, west3$MASE,
  "tuned sarima", west4$MAE, west4$MASE,
  "auto arima", west5$MAE, west5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()
model mae mase
initial arima 62.19197 0.7886397
tuned arima 46.88033 0.5944769
initial sarima 62.19197 0.7886397
tuned sarima 55.78604 0.7074078
auto arima 45.78169 0.5805454